% IRFSIM.m
% 
% Computes the IRF with an announcement of a future change (during the response period) using the model in MODEL.m
%
%                        1       2     3     4         5              6          7       8        9       10 
% Variables in Y(t) = [ yhat(t) pi(t) r(t) ghat(t) E_t(yhat(t+1)) E_t(pi(t+1)) ahat(t) ehat(t) ln(z(t)) zhat(t)]'
% 
%
%                           1      2    3     4
% Variables in eps(t) = [ ea(t) ee(t) ez(t) er(t) ]'
%
%
% Variables in eta(t) = [etay(t) etapi(t)]' 
%
% Date Created: 5-May-2008
% Last Modified: 7-May-2008
%

% clear all
clc
close all


%% Parameters of the IRF

param

n1 = 6;
n2 = 2;
nleads = 1;
s = nleads;

nperiods = T;

if  or(ta < 1, up_to < T+ta) 
    disp('Possible Error: increase time of annoucement, ta')
    disp('Possible Error: increase time of simulation, up_to')
    return
end

%% Load Dimensions
nvars = 10;
nshocks = 4;
nfcast = 2;
%% Load Initial Matrices
% Note: row indicates equation number and column indicates variables number
GAM0 = zeros(nvars);
GAM0(1,1) = 1; 
GAM0(1,3) = 1/sigma;
GAM0(1,5) = -1;
GAM0(1,6) = -1/sigma; 
GAM0(1,7) = -(1-rho_a)/sigma;
GAM0(2,1) = -psi*sigma;
GAM0(2,2) = 1+beta*alpha;
GAM0(2,6) = -beta;
GAM0(2,8) = 1;
GAM0(2,10) = psi;
GAM0(3,1) = -rho_y;
GAM0(3,2) = -rho_pi;
GAM0(3,3) = 1;
GAM0(3,4) = -rho_g;
GAM0(4,7) = 1;
GAM0(5,8) = 1;
GAM0(6,9) = 1;
GAM0(7,1) = -1;
GAM0(7,4) = 1;
GAM0(8,9) = -1;
GAM0(8,10) = 1;
GAM0(9,1) = 1;
GAM0(10,2) = 1;

GAM1 = zeros(nvars);
GAM1(2,2) = alpha;
GAM1(3,3) = rho_r;
GAM1(4,7) = rho_a;
GAM1(5,8) = rho_e;
GAM1(6,9) = rho_z;
GAM1(7,1) = -1;
GAM1(9,5) = 1;
GAM1(10,6) = 1;

C = zeros(nvars,1);
C(1,1) = -(1/sigma)*log(beta);
C(2,1) = (1+beta*alpha-alpha-beta)*pitarget;
C(3,1) = (1-rho_r)*(pitarget-log(beta)) - rho_pi*pitarget;
C(6,1) = (1-rho_z)*log(z);
C(8,1) = -log(z);

PSI = zeros(nvars,nshocks);
PSI(3,4) = 1;
PSI(4,1) = 1;
PSI(5,2) = 1;
PSI(6,3) = 1;

PPI = zeros(nvars,nfcast);
PPI(9,1) = 1;
PPI(10,2) = 1;

%% Load Final Matrices
GAM0n = zeros(nvars);
GAM0n(1,1) = 1; 
GAM0n(1,3) = 1/sigma_n;
GAM0n(1,5) = -1;
GAM0n(1,6) = -1/sigma_n; 
GAM0n(1,7) = -(1-rho_a_n)/sigma_n;
GAM0n(2,1) = -psi_n*sigma_n;
GAM0n(2,2) = 1+beta_n*alpha_n;
GAM0n(2,6) = -beta_n;
GAM0n(2,8) = 1;
GAM0n(2,10) = psi_n;
GAM0n(3,1) = -rho_y_n;
GAM0n(3,2) = -rho_pi_n;
GAM0n(3,3) = 1;
GAM0n(3,4) = -rho_g_n;
GAM0n(4,7) = 1;
GAM0n(5,8) = 1;
GAM0n(6,9) = 1;
GAM0n(7,1) = -1;
GAM0n(7,4) = 1;
GAM0n(8,9) = -1;
GAM0n(8,10) = 1;
GAM0n(9,1) = 1;
GAM0n(10,2) = 1;

GAM1n = zeros(nvars);
GAM1n(2,2) = alpha_n;
GAM1n(3,3) = rho_r_n;
GAM1n(4,7) = rho_a_n;
GAM1n(5,8) = rho_e_n;
GAM1n(6,9) = rho_z_n;
GAM1n(7,1) = -1;
GAM1n(9,5) = 1;
GAM1n(10,6) = 1;

Cn = zeros(nvars,1);
Cn(1,1) = -(1/sigma_n)*log(beta_n);
Cn(2,1) = (1+beta_n*alpha_n-alpha_n-beta_n)*pitarget_n;
Cn(3,1) = (1-rho_r_n)*(pitarget_n-log(beta_n)) - rho_pi_n*pitarget_n;
Cn(6,1) = (1-rho_z_n)*log(z_n);
Cn(8,1) = -log(z_n);

PSIn = zeros(nvars,nshocks);
PSIn(3,4) = 1;
PSIn(4,1) = 1;
PSIn(5,2) = 1;
PSIn(6,3) = 1;

PPIn = zeros(nvars,nfcast);
PPIn(9,1) = 1;
PPIn(10,2) = 1;

[n k] = size(PPI);
[n l] = size(PSI);
%% Initial condition
% Start from the steady state of the original system
[S0_i, S1_i, S2_i, S3_i, unique_i] = Smats(C, GAM0, GAM1, PSI, PPI);
y0 = inv(eye(n) - S1_i)*S0_i;

%% Final steady state
[S0_f, S1_f, S2_f, S3_f, unique_f] = Smats(Cn, GAM0n, GAM1n, PSIn, PPIn);
y0_new = inv(eye(n) - S1_f)*S0_f;

%% Computing IRF to Shock 1--ea(t) 
%  Size of shock is one standard deviation
e1 = zeros(l,up_to);
e1(1,1) = Omega(1,1);

yirf_1 = zeros(n,up_to);

% Store endogenous variables in the transition period
yirf_1_transition = zeros(n,T);

% IRF with initial structure up to period ta-1
yirf_1(:,1) = y0 + S2_i*e1(:,1);
for t = 2:ta-1
    yirf_1(:,t) = S0_i + S1_i * yirf_1(:,t-1);
end

%% LREAnt.m
% Construct sequences of time-varying matrices 
GAM0s = zeros(n*(T+1),n);
GAM1s = zeros(n*(T+1),n);
Cs = zeros(n*(T+1),1);
PSIs = zeros(n*(T+1),l);
PPIs = zeros(n*(T+1),k);

GAM0s(1:n*T,:) = repmat(GAM0,T,1);
GAM0s(n*T+1:end,:) = GAM0n;

GAM1s(1:n*T,:) = repmat(GAM1,T,1);
GAM1s(n*T+1:end,:) = GAM1n;

Cs(1:n*T,:) = repmat(C,T,1);
Cs(n*T+1:end,:) = Cn;

PSIs(1:n*T,:) = repmat(PSI,T,1);
PSIs(n*T+1:end,:) = PSIn;

PPIs(1:n*T,:) = repmat(PPI,T,1);
PPIs(n*T+1:end,:) = PPI;

complex_tol = 1e-6;
% Load anticipated shocks
eps_a = zeros(l,T); 
% Load unanticipated shocks for stochastic simulations.
et = zeros(l,2*T); % size needs to larger than T

% For initial period of the transition
[yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,yirf_1(:,ta-1),eps_a, et(:,ta),complex_tol);
y0 = yseq(:,1);
yirf_1_transition(:,1) = yseq(:,1);
% Loop is useful for stochastic simulation.
for t = 2:T
    eps_a = zeros(l,T-t+1);
    GAM0s(1:n,:) = [];
    GAM1s(1:n,:) = [];
    Cs(1:n,:) = [];
    PSIs(1:n,:) = [];
    PPIs(1:n,:) = [];
   [yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T-t+1,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, et(:,ta+t-1),complex_tol);
    y0 = yseq(:,1);
    yirf_1_transition(:,t) = yseq(:,1);
end

yirf_1(:,ta:T+ta-1) = yirf_1_transition;

%% End of Response with final structure

for t = T+ta:up_to
     yirf_1(:,t) = S0_f + S1_f * yirf_1(:,t-1);
end


%% Transformation of Units
% For shock 1
output_resp_1 = yirf_1(1,:);
g_resp_1 = yirf_1(4,:);
g_resp_1 = 400*g_resp_1;
inflation_resp_1 = yirf_1(2,:)*400;
r_resp_1 = yirf_1(3,:)*400;
inflationE_resp_1 = yirf_1(6,:)*400;
realr_resp_1 = r_resp_1 - inflationE_resp_1;

%% Plot Options 
time = 1:up_to;


